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Time-dependent models of the structure and stability of 
self-gravitating protoplanetary discs. 
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ABSTRACT 

Angular momentum transport within young massive protoplanetary discs may be dom- 
inated by self-gravity at radii where the disk is too weakly ionized to allow the devel- 
opment of the magneto-rotational instability. We use time-dependent one-dimensional 
disc models, based on a local cooling time calculation of the efficiency of transport, 
to study the radial structure and stability (against fragmentation) of protoplanetary 
discs in which self-gravity is the sole transport mechanism. We find that self-gravitating 
discs rapidly attain a quasi-steady state in which the surface density in the inner disc 
is high and the strength of turbulence very low [a ^ 10^'^ or less inside 5 an). Tem- 
peratures high enough to form crystalline silicates may extend out to several au at 
early times within these discs. None of our discs spontaneously develop regions that 
would be unambiguously unstable to fragmentation into substellar objects, though 
the outer regions (beyond 20 au) of the most massive discs are close enough to the 
threshold that fragmentation cannot be ruled out. We discuss how the mass accretion 
rates through such discs may vary with disc mass and with mass of the central star, 
and note that a determination of the M-M* relation for very young systems may allow 
a test of the model. 

Key words: stars: formation — stars: pre-main-sequence — circumstellar matter — 
planetary systems: protoplanetary discs — planetary systems: formation 



1 INTRODUCTION 

Low mass stars fo rm from the collapse of cold, d ense molec- 
ular cloud cores (jTerebev. Shu fc CassenI Il984l ). Although 
the rotation rates of such cores are generally quite small 
ICaselli et al.|[20o3 ) they nonetheless contain amounts of an- 
gular momentum far in excess of the rotational angular mo- 
mentum of a single star. Most of the mass must therefore 
pass through a protostellar disc, and the answers to many 
open problems in star and planet formation hinge on the 
nature of the angular momentum transport that is needed 
for disc accretion. 

In most astrophysical discs the fundamental question 
of what mechanism dominates angular momentum trans- 
port is widely considered to have been solved - MHD turbu- 
lence initiated by the magneto-rota tional instability (MRI) 
can provide the necessary v iscosity l|Balbus fc Hawlevlll99ll : 
iPapaloizou fc NelsonI 120031 '). Protoplanetary disks, however, 
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are so cold and dense that thermal processes proba- 
bly fail to yield even the very small degre e of ion iza- 
tion needed to sustain MHD turbulence ((Blaes fc BalbusI 
1 1994 ). Under these conditions disc self-gravity may pro- 
vide an alternate and possibly dominant mechanism for 
transporting angular mo mentum through the growth of the 
gravitational instability llToonir^ll964l : iLiii fc Pringlelll987l : 



iLaughlin fc Bodenheimedll99i r 



Study of the development of gravitational instability 
in protostellar discs has often been motivated largely by 
the possibility that the instability will lead to fragmen- 
t at io n of the d isc and the formation of gas giant plan- 
ets (|BossI 1 19981 ). The conditions for fragmentation are, 
however, quite difficult to achieve, especially in the in- 
ner, p l anet-forming regions (iMatzner fc Levinll2005l: RafikoyI 
2005; Bolcv ct al. 2006; 'Whitworth fc StamatellosI l2006l : 
StamatcUos fc Whitworth 2008 ; Forgan ct al, 2009). What 
seems more likely is that discs will evolve towards quasi- 
steady states in which the i nstability acts to transport an- 
gular momentum outwards llGammiell2()0l]: iRice et al. I l2003l : 
iLodato fc Ric3l2004 IVorobvov fc Basull2007l ). This is of in- 
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terest in its own right, as it imphes that the conditions 
within young protoplanetary discs - at precisely the epoch 
when planetesimals and perhaps larger bodies are forming 
- may largely be set by the physics of angular momentum 
transport via gravitational instability. The recent progress in 
understanding the conditions for fragmentation has yielded 
a much clearer understanding of this physics, which we use 
here to construct realistic time-dependent models for the 
structure of self-gravitating protoplanetary discs. Our mod- 
els rely on two specific properties of transport via gravita- 
tional instability. First, that transport can be approximated 
as a l ocal viscous process for al l except the most massive 
discs l|Lodato fc Ricell2004 | 200^ . Second, under conditions 
of thermal equilibrium the strength of angu lar momentum 
transport is set b y the cooling rate of the disc ()Gammiell200ll : 
iRice et ahlbOOST ). Our results show that self-gravitating discs 
will settle into quasi-steady states that appear independent 
of the initial conditions, but that their properties are quite 
unlike those that result from angular momentum transport 
via generic turbulent processes. Specifically, we find that the 
surface density profile is reasonably steep and in the cases 
considered here, ~ 80 % of the mass within 50 au is lo- 
cated inside 10 — 20 au. The quasi-steady mass accretion 
rate depends strongly on the disc mass and on the mass 
of the central star. For a constant star to disc mass ratio, 
however, the relationship between mass accretion rate and 
central star mass is similar to that found observationally.We 
further show that the secular evolution of the disc does not 
typically result in any regions that would be unstable to 
fragmentation, with the region inside 10 — 20 au being par- 
ticularly st able unless some mechanism - such as convection 
(|Bosj|2004l ) - can significantly increase the cooling rate. 

Our focus in this paper is on the outer cool regions 
of protoplanetary discs, and for this reason and for sim- 
plicity we consider disk models in which self-gravity pro- 
vides the sole source of angular momentum transport. Of 
course both the extreme inner region (inside 0.1 au) and 
the upper layers of the disc further out could become suffi- 
ciently ionised for the MRI to operate (jGammidliggel ). and 
this would modify our results and admit new physical ef- 
fects. In particular, a pile-up of material brought to within 
1 — 2 au by self-gravity could eve ntually trigger the onset of 
the MRI and episodic outbursts (lArmitage. Livio fc Pringld 



2.1 Viscous evolution 



l200ll : IZhu. Hartmann fc Gammiell2009 '). These might be re- 
lated to the FU Orionis phenomenon (Hartmann fc Ken von 
1 19961 ). Ultimately the evolution of discs around very young 
stars could be driven by both MRI and the gravitational 
instability (|Terauemll20CI8l '). 

The paper is organised as follows. In Section 2 we de- 
scribe how we can self-consistently model protostellar discs 
evolving through self-gravity alone. In Section 3 we describe 
our results, and in Section 4 we summarise our conclusions. 



2 DISC MODELS 

We model the evolution of a self-gravitating protoplanetary 
disc under the assumptions that the potential is fixed, that 
the disc is in thermal equilibrium at all radii, and that the 
local strength of angular momentum transport is set by self- 
gravity. 



The surface density S(r, t) of an axisymmetric disc evol ves 
according to (|Lvnden- Bell fc Pringi3ll974l : |Pringldll98lD 
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where v is the kinematic viscosity. If the disc is able to attain 
a steady-state this equation can be integrated to give an 
expression for the mass accretion rate M which, a t radii 
large compared to the radius of the star, is (IPringldUgSl. ) 



M = 37r!/E. 



(2) 



The steady-state radial surface density profile is therefore 
determined by t he radial profile of the vis cosity. We adopt 
the a formalism l|Shakura fc Sunvaevll973l ) in which the vis- 
cosity is taken to be = aCsH where Ca is the disc sound 
speed, H = Cs/fl is the disc scale height, and a << 1 is a 
parameter that determines the efficiency of angular momen- 
tum transport. 

The strength of self-gravity depends upon the thermal 
properties of the disc. In thermal equilibrium the viscosity 
generates dissipation i n the disc a t a rate D{R) per unit area 
per unit time, where (iBell fc Lin|[l994 ) 



D{R) 



(3) 



which must balance cooling via thermal emission from the 
disc surfaces (if the disc is optically thick) or via optically 
thin emission integrated vertically through the disc. 



2.2 Determination of a 

An accretion disc can become gravitationally unstable if 
(|Toomrdll963 l 



ttGE 



(4) 



where k is the epicyclic frequency which is replaced by the 
angular frequency, 57, in a Keplerian disc. One possible out- 
come is that unstable discs fragment to produce bound ob- 
jects and has been su ggested as a poss ible mechanism for 
forming giant planets (|BossI [l998l . |2003 ). For axisymmetric 
instabilities this requires Q < 1, while for non-axisymmetric 
i nstabilities this can occur for Q values as high as 1.5 — 1.7 
(|Durisen et al.|[2007l ). It has, however, been realised recently 
that the above condition is not sufficient to guarantee frag- 
mentation. The disc must also be able to l o se the energy gen- 
erated by the instability l|Gammid [2OO1I : iRice e t al. 2003). 
The rate at which energy mus t be lost depends on the equa- 
tion of state l|Rice et al. l2005l ). but in protostellar discs, the 
cooling time w ould generally need to satisfy tcooi ^ 3f2~^ 
(|Gammidl2"00ll ). 

For cooling times greater than that required for frag- 
mentation, the disc will settle into a quasi-steady state in 
which the instability acts to transport angular momentum 
llLaughlin fc Rocvzkal 1 19961 : iGammid I2OOII : iLodato fc Rlcd 
1200 ^. In principle this transport need not be local 
(iBalbus fc Papaloizoulll999l '). and if this were the case then 
neither writing a as a function of local conditions nor using 
equation [T] would be valid. Simulations, however, show that 
the local approximation is surprisingly good at capturing 
the behavior of self-gravitating protoplanetary discs, which 
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can therefore be regarded as having an effective viscosity of 
the form v = acsH. Moreover, under conditions of thermal 
equihbrium the value of a can be derived via a simple energy 

balance requirement. 

The disc cools at a rate (jPringlel Il98ll : 
IJohnson fc Gammiell2003l ) 



A = 2aT^ 



(5) 



where a is the Stefan-Boltzmann constant, Te is the effective 
temperature, and the factor of 2 comes from the radiation 
escaping on both sides of the disc. If the disc is optically 
thick, radiation transport can be treated in the diffusion 
approximation, and the effect ive temperature can be shown 
to be given by l|Hubenvlll990l ) 



8Tl_ 
3 r ' 



(6) 



where Tc is the temperature of the disc midplane, and r is 
the Rosseland mean optical depth. The internal energy per 
unit area, U, is then 



U 



7(7 ■ 



(7) 



where 7 is the specific heat ratio. 

The coohng time is then simply tcooi ~ U/A, where 
A is the cooling function given by equation ([5}. Since we 
are assuming that the disc is in a quasi-steady state, the 
cooling must be balanced by dissipation. If we assume that 
the energy is dissipated locally, w hich app ears to be the cas e 
for Q ~ 1 (jBalbus fc Papaloizoul ll999: Lo dato fc Ricel2004h . 
we can use equation ^ to get |Gammiell200ll ) 



(8) 
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We should acknowledge that the nature of the transport - 
whe ther loca l or global - may depend on the form of the cool- 
ing (jMeiia et a l. 2005: ,Durisen et al.|[2007il and that our use 
of the local approximation is clearly a simplication. There 
may well be situations in which the energy is not dissipated 
locally and although the disc as a whole may be in thermal 
equilibrium, equation (|8]) is not strictly valid. 



2.3 Modelling a fully self-gravitating disc 

To consider the evolution of a disc in which angular mo- 
mentum transport is governed by self-gravity alone, we as- 
sume that the disc settles into a quasi-steady state and that 
energy is dissipated locally. We assume that in this quasi- 
steady state the Toomre Q parameter satisfies Q = 1.5, un- 
less Tc < 10 K in which case we set Tc = lOK. For a given 
surface density, equation Q can be used to determine the 
midplane sound speed, Cs- The midplane temperature can 
then be determined using 



Tc 



"fkB 



(9) 



where fi = 2.4 is the mean molecular weight, irip is the 
proton mass, and fcs is Boltzmann's constant. 

For a given surface density, E, and using the assumption 
that Q = 1.5 (unless T < lOK) we can now determine the 
midplane temperature. To determine the cooling rate, how- 
ever, we need to determine the effective temperature, Te, for 
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Figure 1. Opacities from lBell fc LirJ lll994l') against temperature 
for densities of p = IQ-^, 10"*, 10"'', IQ-i^, IQ-^, 10"", IQ-^ g 



which we need the optical depth. We can approximate the 
optical depth using 



r = I dzK.{pz,Tz)pz ~ Hk{p,T)p, 
lo 



(10) 



where k is the Rosseland mean opacity and p and T are 
an average density and temperature for which we use p = 
E/(2iif) and T = Tc- For the Rossela nd mean opacitie s 
we use the analytic approximations from lBell fc LinI l|l994 ). 
Figure [1] shows the opacity against temperature for a range 
of densities. The opacity is dominated, in order of increas- 
ing temperature, by ice grains, metal grains. Molecules, H~ 
scattering. Bound-free and free-free absorption, and electron 
scattering. The sudden changes in opacity at T ~ 170 K and 
T ~ 1000 K are due to the evaporation of the ice mantles 
on the ices grains and the evaporation of the metal grains 
respectively. 

Once the optical depth, r, is known the cooling func- 
tion, A, is calculated using a modified form of equations ([6} 
and ©, 



T„*) 



1 +t2 



(11) 



where To = lOK is assumed to come from some background 
source that prevents the midplane cooling below this value 
(jStamatellos et al]l2007bl l , and the last term is introduced to 
smoothly int erpolate between optically thick and optically 
thin regions (| Johnson fc Gammid [2OO3I ) . The cooling time 
is then tcooi = U / A. with U given by equation ([7} and the 
viscous a can then be determined using equation ([8]). The 
viscosity is then u = otCsH and the disc is evolved using 
equation ([TJ. 



3 RESULTS 

3.1 Surface density profiles 

The simulations presented here assume, initially, a rather 
extreme case of a central star with a mass of 0.35 M© sur- 
round by a circumstellar disc with a mass of 0.35 M©. The 
disc is assumed to have an initial power-law surface density 
profile E oc between 0.5 and 50 an. The assumption that 
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1 10 
R (au) 

Figure 2. Figure showing how 0.35 Mq discs around a 0.35 Mq 
star with initial surface density profiles of S oc (top panel), 
S oc r~^'^ (middle panel), and S oc (bottom panel) evolve 
into quasi-steady states. In all the cases the resulting quasi-steady 
profiles are the same. 



Q = 1.5 then determines the disc midplane temperature, Tc, 
with the additional constraint that Tc ^ 10 K. The process 
described above determines a and the disc is then evolved 
in time using equation ([T}. We consider 3 different initial 
surface density profile, /3 = 2.0, /3 — 1.5, and /3 = 1.0. 

The surface density evolution for the three 3 initial pro- 
files is shown in Figure [2] We should stress that we do not 
assume that the initial profiles are necessarily realistic, but 
simply want to if establish if such discs will settle into quasi- 
steady states and if such a state depends in any way on 
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Figure 3. Figure showing the quasi-steady surface density profile 
for a 0.35 Mq disc around a 0.35 Mq star (solid line). The profile 
moves inwards with time (dashed line) as mass accretes onto the 
central star. 



the initial surface density profiles. Each panel in Figure [5] 
shows the initial power-law surface density, 4 subsequent 
surface density profiles each separated by 20000 years and 
the final surface density (dashed-line) at the end of the sim- 
ulation which we stop after 10® years. What this shows is 
that the surface density adjusts in all three cases to a state 
that appears to be largely independent of the initial sur- 
face density profile. For (3 — 1.0 and l3 = 1.5 this settling 
takes a some time (~ 80000 years), while for /3 = 2, it oc- 
curs almost immediately. For clarity, Figure|3]shows the sur- 
face density profile when a quasi-steady state has just been 
achieved (solid line). With time, this profile moves inwards 
and the disc mass decreases, as illustrated by the dashed 
line in Figure[3]which shows the surface density profile when 
Mdisc = 0.25 Mq. 

In all cases what is happening is that mass is redis- 
tributing itself to produce a state in which the accretion 
rate, M, is largely independent of r. The is illustrated in 
Figure |3] which shows the evolution of the mass accretion 
rate for /3 — 1.5. The behaviour is essentially the same for 
/3 = 1 and {3 = 2 except that the latter reaches an approxi- 
mately constant mass accretion rate more rapidly than the 
other two cases. In the prescription presented here, the vis- 
cous a depends on the local cooling rate - which depends on 
the local temperature - and the accretion rate depends on 
the local viscosity and surface density (e.g., equation (O). 
For initially fiatter surface density profiles (S oc and 
E cx r~^'^) the outer disc contains most of the mass. To 
be gravitationally stable (see equation |4|, the temperature 
in the outer disc also needs to be high. The inner disc, on 
the other hand, can have relatively low temperatures and 
remain gravitationally stable. The cooling time, and conse- 
quently viscosity, in the outer disc is therefore significantly 
greater than in the inner disc, producing an initially much 
larger mass accretion rate in the outer disc than in the in- 
ner disc. This process, however, moves mass inwards and in 
doing so not only increases the surface density in the inner 
disc, but also the temperature to keep the disc gravitation- 
ally stable. This increases the cooling rate in the inner disc, 
the viscosity, and the mass accretion rate. Eventually the 
mass is redistributed such that the entire disc has the same 
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mass accretion rate. If the initial surface density profile is, 
however, initially close to the quasi-steady profile - as in the 
third panel of Figure[2]- very little redistribution is required 
and the disc very quickly settles into a quasi-steady state. 

This mass redistribution leads to quasi-steady sur- 
face density profiles that are reasonably steep (E ^ '"~^) 
and that have quite a substantial break at 10 - 20 au. 
Even though the discs extend to 50 au, most of the mass 
(> 80 %) is actually located in the inner disc. It is now 
generally a ccepted that discs disappear within a few mil- 
lion years (|Haisch. Lada fc Ladal I2OOII) and consequently 
that gas giant planets need to form within this timescale. 
A substantial reservoir of mass in the inner disc could 
significantly accelerate planet growth if ga seous planets 
form via the standard core accretion scenario l|Pollack et al.l 
1 19961 ). This type of profile is also qualitatively consistent 
with observations of massive discs ([Rodriguez ct al. 2005; 
ICarrasco-Gonzalez et al]|2009l ) that suggest that such discs 
have small radii (rdisc ~ 25 au). It is also quite likely that 
with so much mass located in the optically thick inner re- 
gions, and the sub-mm fiux largely determined by the cold, 
outer disc, the masses of such di scs could be easily under- 
estimated (|Hartmann et al]|2006l ). Such gas rich inner disc 
could be detected using near-IR gas emission lines such as 
the CO bandhead (|Thi et al.ll2005l : lNaiita et al.ll2007l '). How- 
ever, for the disc masses considered here, these lines are still 
optically thick and will still not provide accurate probes of 
the disc mass. 

Although the goal here is to illustrate the quasi-steady 
structure of a self-gravitating disc, it is also interesting to 
consider the time evolution of these systems. Firstly, since 
the quasi-steady surface density profile is independent of the 
initial conditions, continuing these simulations to lO'' years 
allows us to compare the quasi-steady nature for different 
disc masses. There is also no reason to suspect that the mass 
falling onto the disc will do so in such a way as to immedi- 
ately produce a quasi-steady disc. We might therefore expect 
that the disc will not initially be in equilibrium. The settling 
times achieved here (< 10^ years) are similar to an d proba- 
bly le ss than typical free-fall times (Ward- Thomps on et al.l 
I2OO7I ). suggesting that these systems could quite quickly at- 
tain quasi-steady states. Figure [2] also shows that beyond 
~ 1 au the disc reaches a quasi-steady state in 20000 years 
or less and might imply that these systems are rarely out of 
equilibrium. 

3.2 Temperature profiles 

The evolution of the midplane temperature is shown in Fig- 
ureO Again, the Figure shows the initial temperature profile 
followed by 4 subsequent temperature profiles separated by 
20000 years, and a final temperature profile (dashed-line) 
after 10® years. Since the temperature is specified using 
Q = 1.5 and Tc ^ 10 K, and since the quasi-steady pro- 
files are essentially the same in all cases, we show only the 
temperature profile for E oc r~^'^ . From Figure[5]and Figure 
[T]we can start to understand the disc structure. The outer 
disc has Tc < 170 K and so the opacity is dominated by ice 
grains. The asymptotic nature of the temperature profile at 
large radii is simply due to the assumption that Tc ^ 10 K. 
The change in slope at, initially, a radius between 10 and 
20 au that moves inwards with time is due to the evapora- 
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Figure 4. Figure showing the mass accretion rate against radius 
at different times for S oc r~^-^. The almost diagonal line shows 
the initial mass accretion rate while the other solid lines show the 
accretion rate at 20000 year intervals until a quasi-steady rate is 
achieved with a mass accretion rate of M ~ 10~^ Mq/jt. The 
dashed hne shows the accretion rate after 10® years when the 
disc mass has decreased to 0.25 M©. The two other initial surface 
density profiles (which we don't show here) evolve to the same 
quasi-steady accretion rate, but over slightly different timescales, 
with the E oc r~'^ disc evolving to a quasi-steady state much more 
rapidly than S oc and E oc r~^'^ discs. 
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Figure 5. Disc Midplane temperature against radius for E oc 
j,-i.5 rpj^g temperature profiles for E oc and E oc aren't 
shown since the quasi-steady profiles are the same. The figure 
shows the initial temperature profile followed by 4 profiles sepa- 
rated by 20000 years, and the profile after 10® years (dashed line). 
The form of the temperature profile is essentially determined by 
the opacity. Beyond 10 au, ice grains dominate the opacity. These 
melt at 170 K, producing the change in profile at ~ 10 au. Within 
10 au, metal grains dominate the opacity until the temperature 
reaches 1000 — 2000 K at which point these grains evaporate. 



tion of the ice mantles on the grains and is essentially the 
"snowline" . The opacity of the inner part of the disc (within 
~ 10 au) is dominated by metal grains until the tempera- 
ture reaches 1000 — 2000 K where the grains evaporate and 
the temperature cannot rise significantly due to the sudden 
decrease in opacity. This results in a further change in the 
surface density profile. The quasi-steady nature is therefore 
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quite strongly dependent on the opacity of the disc material. 
In the vertically averaged approximation the optical depth 
is evaluated using only the midplane value of the opacity. 
In the hot inner regions where the midplane temperature 
exceeds the dust sublimation temperature this is likely to 
result in an underestimation of the true optical depth, since 
the surface regions are cooler and could have a larger opac- 
ity. To evaluate the true optical depth in detail would require 
solving for the vertical structure and considering the time 
scales for the evaporation and condensation of particulates 
within the turbulent flow, which we do not attempt. Clearly, 
however, the sense of the effect will be to increase the opti- 
cal depth and cooling time still further and hence the disc 
temperature at the midplane will remain above the grain 
evaporation temperature. 

As will be discussed in more detail later, the a values 
in the inner regions of these discs are low and, as shown in 
Figure [S] within ~ 3 au are < 10~*. It is often assumed that 
a ~ 10~^ which, for mass accretion rates similar to those 
shown here, would imply much lower temperatures in the in- 
ner disc. Crystalline silicates - w hich are commonly observed 
in pr otostellar discs ( Bouwman e t al. 2001; van Bockcl ct ahl 
|2004 ) and require T > 800 K - are therefore generally 
thought to for m close to the central star (r ~ 0.1 au for 
TTauri stars) l|Dullemond. Apai fc WalchI |2006| ) and then 
transp orted outw ards via turbulent mixing to at least 10 - 
20 au (|Gailll200ll ). Figure[Sl however, suggests that when the 
transport processes are dominated by self-gravity, the tem- 
perature could be high enough to form crystalline silicates 
at radii of a few au. 



3.3 Mass Accretion 

As discussed above (and illustrated in Figure the disc 
relatively quickly settles into a quasi-steady state with an 
approximately constant mass transfer rate. The correspond- 
ing a values are shown in Figure [S] Again this Figure shows 
the initial a profile, 4 subsequent a profiles separated by 
20000 years, and the final a profile after 10^ years. The 
strong dependence of a on the midplane temperature means 
that small variations in Tc can produce large variations in a. 
This only really occurs in the outer disc when Tc is close to 
the minimum temperature of 10 K, or where opacity changes 
produce sudden changes in Tc coupled with the initial condi- 
tions being far from the equilibrium conditions. The a and 
M profiles have therefore been averaged to remove these 
small-scale variations. What Figure [6] illustrates is that the 
a value in a quasi-steady, self-gravitating disc is not con- 
stant and can reach extremely low-values in the inner disc. 
In a more complete model other sources of transport (such 
as weak turbulence in the disc mid-plane excited by MRI 
active zones at the disc surface) might well dominate in this 
region. Such a process could also result in episodic accretion 
events if the other source of transport were to depend on the 
disc properties - such as temperature in the case of MRI - 
and could explain phenomen a such as FU Orionis outbursts 
l|Hartmann fc KenvonI [l996l V The mass accretion rate into 
the inner disc and ultimately onto the star, however, would 
still be set by the self-gravitating t ransport taking place at 
larger radii (|Bolev fc Durisenll2008l ). 

Figure |4] also shows that for our chosen parameters 
(Af* = 0.35 Mq and Mdisc = 0.35 M©) the initial quasi- 
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Figure 6. Viscous a against radius at different times for S oc 
J, -1.5 'j'jjg results for S oc and S oc are essentially 

the same, so aren't shown. The figure shows the initial a profile 
followed by 4 a profiles separated by 20000 years (solid lines) until 
a quasi-steady state, and the a profile after 10^ years (dashed 
line). It has been shown that fragmentation can occur if a ~ 0.06. 
The a values in the inner disc can be very small, well below that 
required for fragmentation. The a value in the outer disc can, 
however, become quite large (> lO"'^) but in this case is still 
below - but close to - that required for fragmentation. 



steady mass accretion rate is ~ 2 x 10"'' M0/yr. The dashed 
line in Figure |4] shows the accretion rate after 10® years. At 
this stage, the disc mass is ~ 0.25 Mq and the accretion 
rate is ~ 10~* Mq/jt. Since this is the quasi-steady rate for 
a 0.25 Mq disc around a 0.35 Mq star, this illustrate that 
the accretion rate depends non-linearly on the disc mass (a 
factor of 1.4 reduction in disc mass reduces the accretion 
rate by about an order of magnitude). 

This dependence of the quasi-steady accretion rate on 
the disc mass can be roughly understood by considering the 
method discussed in section 2. For a constant Toomre sta- 
bility parameter, Q, and central mass, the sound speed, Cs, 
depends linearly on the surface density. Since Cs oc x/TT, this 
gives Tc oc E'^. Using u = aCs/f2, the mass accretion rate - 
shown in equation ([2j - can be rewritten as AI oc aE^. Com- 
bining equations ((7}, (|8]), (|10p . and HI}, and using r « kE, 
gives a oc E* and consequently oc E''. Here, for simplic- 
ity, we ignore To in equation pi|l and assume that the disc 
is optically thick so r/(l -|-r^) « r. As E is essentially deter- 
mined by the disc mass, the accretion rate then depends on 
the disc mass to the seventh power. A factor of 1.4 reduction 
in disc mass would produce an order of magnitude change 
in accretion rate, consistent with that shown in Figure ID 



3.4 Fragmentation 

It has been shown (iGammid [2OO1I : iRice et al.1 [200I that 
there is a minimum cooling time, tcooi, for which a self- 
gravitating disc can remain in a quasi-steady state with- 
out fragmenting to form bound objects. The exact value 
of the minimum cooling time depends on the equation 
of state (|Rice et al.ll2005f) with fragmentation occ uring for 
tcooi ^ 3il~^ when the s pecific heat ratio 7 = 5/3 (jCammid 
|200l| ). lRice et al.1 (|2005D . however, show that for all values 
of 7 the fragmentation boundary occurs for a ~ 0.06. Fig- 
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ure [6] shows that once the disc settles into a quasi-steady 
state, the a values in the inner disc (r < 10 au) are well 
below that required for fragmentation. Between 20 and 50 
au, however, a is almost constant and has a value {a ~ 0.04) 
just below that required for fragmentation. A small increase 
in disc mass would, however, increase a and could produce 
conditions that wo uld allow for fragmentation in the outer 
parts of such discs (jStamatellos et al]|2007al : iGreaves et al.l 
|2008|). We should, how ever, add t h at th e results obtained 
bv lGammiel l|2001f ) and iRice et al.l (|2003l ) were for systems 
in wh ich the cooli ng time was chosen to be the same at all 
radii. [Meiia et ahl (2005) argue that in systems with a radi- 
ally varying cooling time, the fragmentation boundary may 
not be a simple function of cooling time and that there could 
be regions with a > 0.06 that do not fragment. 



3.5 Variation with central mass 

We have also considered how the mass of the central object 
influences the quasi-steady nature of a self-gravitating disc. 
To do this we consider a disc mass of 0.25 Mq around central 
objects with masses of M, = 0.25 Mq, M* = 0.35 M©, and 
M* = 0.5 Mq. In all cases we assume an initial surface den- 
sity profile of E oc r^^ and evolve the system until a quasi- 
steady state is reached. The resulting quasi-steady surface 
density profiles are shown in Figure [7] with the solid line for 
M, = 0.25 Mq, the dashed line for M, = 0.35 Mq, and the 
dashed-dot line for M» = 0.5 Mq. The surface density pro- 
files are, as expected, quite similar. The slight differences are 
a consequences of the dependence of the Toomre Q parame- 
ter on the mass of the central object. This can be understood 
by considering Figure |8] which shows the temperature profile 
for the three cases. Since we assume Q = 1.5, as the cen- 
tral object mass decreases (decreasing the orbital frequency 
Q.) the sound speed, Cs, and hence temperature, needs to 
increase. This changes the locations where the ice mantles 
melt and where the metal grains evaporate and hence results 
in different central object masses producing slightly different 
- although similar - surface density profiles. 

Figure |9] shows the quasi-steady mass accretion rates 
for A/* = 0.25 Mq (sofid line), M, = 0.35 Mq (dashed line) 
and M, = 0.5 Mq (dash-dot line). There is a very strong 
dependence on the central mass with the accretion rate for 
M* — 0.25 Mq being more than an order of magnitude 
greater than that for M, = 0.5 Mq. Again, we can under- 
stand this dependence by considering the method described 
in section 2. For a constant disc mass, and hence surface den- 
sity (E), a constant Toomre Q requires Cs oc l/fi oc 1/\^M, 
giving Tc oc 1/M,. Using v = aCg/Q, the mass accretion 
rate can be rewritten as M oc aTc''^. Equations ([7]), (|8]), 
(|10p . and (|11|) can again be used to show that a oc Tc^^, 
giving a mass accretion rate of M oc T^^, or M oc l/M* . A 
factor of 2 increase in the mass of the central star therefore 
reduces the mass accretion rate by a factor of 22, consis- 
tent with that seen in Figure|9] If, however, we combine this 
relationship (M oc 1/Mf) with the relationship between ac- 
cretion rate and disc mass (M oc oc Mjig^.), the accretion 
rate actually increase as M oc M« if a constant disc to star 
mass ratio is maintained. This is consistent with the obser- 
vations that in the Orion Nebula clust er massive discs ar e 
more common around lower-mass stars (jEisner et al.ll200^ . 
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Figure 7. Quasi-steady surface density profiles for Mtjigc = 0.25 
Mq and M, = 0.25 Mq (solid line), M* = 0.35 Mq (dashed line), 
and M* = 0.5 Mq (dash-dot line). The profiles are almost the 
same, but with a few differences due to the different temperature 
profiles resulting from the Toomre Q parameter's dependence on 
the mass of the central star. 



Something to point out is that as the disc mass de- 
creases, it become increasingly difficult to sustain steady, 
self-gravitating accretion in the outer regions of the disc. The 
surface density (Figures [2] and [Sj and temperature profiles 
(Figures |S] and O show that most of the mass is in the inner 
disc and that, as the disc mass decreases, the temperature 
in the outer disc approaches 10 K, our assumed minimum 
temperature. Equation (|lip shows that as Tc approaches 10 
K, the cooling function becomes very small, the cooling time 
becomes very long, and the effective a becomes very small. 
Decreasing the disc mass much below 0.25 Mq effectively 
turns off self-gravitating accretion beyond ~ 30 au. This 
is, of course, due to our assumption that the quasi-steady 
Toomre Q value is 1.5. If quasi-steady transport can occur 
for larger values of Q, self-gravitating transport could still 
operate effectively at lower disc masses. This also has im- 
plications for the productio n of brown dwarfs and low -mass 
companions at large radii l|Stamatellos et af] l2007al ) , sug- 
gesting that fragmentation at large radii is unlikely when 
the disc is in a quasi-steady state. This doesn't proclude 
fragmentation at large radii, but implies that it must either 
occur when the disc is extremely massive, or must occur 
when the disc has not yet reached a quasi-steady state. 

The corresponding a profiles are shown in Figure 1101 
The dependence of a on M* discussed above gives larger a 
values for M* = 0.25 Mq than for M* = 0.5 Mq. Again 
the a values in the inner disc are well below those required 
for fragmentation, while those in outer disc are close to, 
but below, the required value (a = 0.06). As discussed ear- 
lier, however, small changes in some of the parameters could 
produce conditions suitable for frag mentation in t he outer 
parts of the disc, but as suggested bv lRafikovl(|2005h . it is ex- 
tremely difficult to see how fragmentation can occur within 
10 au even for the relatively massive discs considered here. 
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Figure 8. Quasi-steady temperature profiles for Mjis;. = 0.25 
M0 and Af, = 0.25 Mq (solid line), M, = 0.35 Mq (dashed 
line), and Af» = 0.5 Mq (dash-dot line). The difference in the 
profiles is a result of the dependence of the Toomre Q on the 
mass of the central star (through Q). As the central star mass 
decreases, the sound speed - and hence temperature - needs to 
increase to maintain a constant Q value. This means that the 
radii at which ice mantles and metals evaporate differs for the 
different central star masses, and produces (as seen in Figure (Tjl 
slightly different surface density profiles. 
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Figure 9. Quasi-steady mass accretion rates for M^i^^ = 0.25 
M0 and M* = 0.25 Mq (solid line), M* = 0.35 Mq (dashed 
line), and Mt = 0.5 Mq (dash-dot line). The accretion rate for 
M* = 0.25 Mq is significantly great than for Mt = 0.5 Mq. This 
is a consequence of the increased disc temperature required to 
keep the disc around the lower mass star gravitationally stable. 
The strong dependence of a on temperature results in the M* = 
0.25 Mq system having an accretion rate more than an order of 
magnitude greater than the M« = 0.5 Mq system. 



4 DISCUSSION AND CONCLUSIONS 

We consider here the evolution of protostellar discs in which 
angular momentum transport is driven only by self-gravity. 
We are particularly interested in the quasi-steady nature of 
such systems, which we define here as a state in which the 
mass accretion rate is approximately the same at all radii. 
We assume that the disc will cool to a state of marginal 
gravitational stability {Q ~ 1.5) and will then be in thermal 
equilibrium. The evolution of such systems is determined by 



10"' 



a 




R (au) 

Figure 10. Effective viscous a profiles for Mjisc = 0.25 Mq 
and Ah = 0.25 Mq (solid line), = 0.35 Mq (dashed line), 
and Alt = 0.5 Mq (dash-dot line). As discussed in the text, for 
a fixed disc mass a increases with decreasing star mass. In all 
cases, however, the a values in the inner disc are well below that 
required for fragmentation. Although larger in the outer disc, the 
values are close to, but still below, that needed for fragmentation. 



the effective gravitational viscosity which is calculated from 
the assumption that the viscosity dissipates energy at the 
same rate as energy is lost from the disc via radative cooling 
- which depends on the chosen op acity. This is in some sense 
the inverse of a " classical" a disc (jShakura Sz Sunvaevll973h 
in which the prescribed a value determines the viscosity, 
which then sets the dissipation and cooling rates. 

We carried out a number of different simulations, ini- 
tially considering a central object of 0.35 Mq surrounded by 
an equal-mass disc with a power-law surface density profile 
from 0.5 au to 50 au. We also considered central masses of 
0.25 Mq and 0.5 Mq with a disc mass of 0.25 Mq. The discs 
are in all cases evolved until a quasi-steady state is reached. 
The primary results are summarised below. 

• The discs reasonably quickly settle into a quasi-steady 
state that is largely independent of the initial surface den- 
sity profile. The settling time does depend somewhat on the 
initial profile, but beyond ~ 1 au a quasi-steady state is 
reached in 20000 years or less, suggesting that such systems 
will rarely be out of equilibrium. 

• The quasi-steady surface density profile is relatively 
steep (S ~ '"~^) with changes in profile that correspond 
to the melting of ice grains (effectively the "snowline") and 
the evaporation of metal grains. The steep profile means 
that most (~ 80 %) of the mass is located within 10 — 20 au, 
even though the discs extend initially to 50 au. Having a lot 
of the mass in the inner disc could aid planet formation, and 
is consistent with curr ent observations showin g that massive 
discs have small radii (|Rodriguez et al]|2005l ). 

• The corresponding temperature profile indicates that 
in a quasi-steady state the inner disc (within ~ 3 au) can 
be hot (Tc > 1000 K) which suggests that such regions 
could be sufficiently ionised to be unstable to the growth of 
MRI (|Balbus fc Hawlevl[l99ll ). The MRI viscosity is hkely 
to be significantly larger than the effective gravitational vis- 
cosity, so will drain the i nner disc, potentially producing 
FU Orionis-like outbursts (jArmitage. Livio fc Pringlg|200ll : 
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IZhu, Hartmann fc Gammiell2009l ). This process could repeat 
once the inner disc has been replenished by material from 
the outer disc. The hot inner disc (Tc > 800 K) also sug- 
gests that the formation of crystalline silicates could occur 
at relatively large radii (r < 3 — 4 au) during the earliest 
stages of star formation. 

• The mass accretion rate in quasi-steady discs also 
has a strong dependence on disc mass and on the mass 
of the central star. From the simulations, and from an- 
alytic approximations, we find M oc oc M^^^^ and 
M oc 1/Mf. Together this suggests that for a constant disc 
to star mass ratio M cc Ml, consi stent with suggestions 
IIKratter. Matzner fc Krumholjr2008l ) that the gravitational 
instability is the most important mechanism for stars with 
final masses > 1 M©. This is also intriguingly similar to 
the M oc Ml relationship between accre tion rate arid ste l- 
lar mass determined observationally ( Muzerolle et al.ll2005h , 
although the simulation results probably apply to systems 
considerably younger than those observed. 

• It has been shown (|Gammiel I2001I : iRice et al.l |2003| ) 
that fragment ation requires rap id cooling and is equivalent 
to a ^ 0.06 l|Rice et al.l l2005h . Despite the discs consid- 
ered here being massive, none of the quasi-steady discs had 
a > 0.06 at any radii. In the outer disc (r > 20 au) the 
a value was close to the fragmentation limit, but in the 
inner disc it was orders of magnitude below that required. 
Small changes in some disc proper ties could lead to fragmen- 
tation conditions a t large radii (IStamatellos et al] l2007al : 
iGreaves et al.|[2008l ). but this does not appear to be pos- 
sible in the inner d i sc, consistent with other calcul ations 
ijRafikovl 12001 l2007l : Iwhitworth fc Stamatellojiiooi l. The 
quasi-steady surface density profile also has most of the mass 
in the inner disc, suggesting that if discs do fragment at large 
radii they either do so prior to settling into a quasi-steady 
state, or do so when the disc is extremely massive. 



Overall, the results presented here suggest a scenario 
in which - while mass is falling from the envelope onto 
the disc - the disc mass remains comparable to the cen- 
tral star mass and the system evolves towards a quasi- 
steady state with mass accretion rates in excess of 10~^ 
M0/yr. MHD turbulence (MRI) will probably also play a 
role i n some reg ions of the disc (|Armitage. Livio fc Pringlg 
I2001I : iTerauemI [2008) and, in particular, could lead 
to FU Orionis outburs ts if the inner disc tempera- 
ture exceeds ~ 1400 K IIArmitage. Livio fc Pringld [2OO1I : 
IZhu. Hartmann fc Gammid \2QQ% 7 Once envelope infall 
ceases, the disc mass will start to decrease, rapidly reducing 
the mass accretion rate due to self-gravity. The midplane 
in the inner disc will have a very low viscosity (a < lO"**) 
and will effectively become dead. This can occur for rela- 
tively large disc masses (Mdisc ~ 0.1 M©) which, together 
with most of the mass being in the inner disc, could enhance 
planet formation. The low viscosities in the midplane of the 
inner disc also mea ns that embedded pla nets can easily clear 
a large inner gap (jSver fc Clarkdll995l l and could explain 
some of the obs erved systems with near-infrared deficits 
l|Hartmannll2008l ) . often referred to as transition systems. 
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